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A practical finite temperature theory is developed for the superfluid regime of a weakly interacting Bose 
gas in an optical lattice with additional harmonic confinement. We derive an extended Bose-Hubbard model 
that is valid for shallow lattices and when excited bands are occupied. Using the Hartree-Fock-Bogoliubov- 
Popov mean-field approach, and applying local density and coarse-grained envelope approximations, we arrive 
at a theory that can be numerically implemented accurately and efficiently. We present results for a three- 
dimensional system, characterizing the importance of the features of the extended Bose-Hubbard model and 
compare against other theoretical results and show an improved agreement with experimental data. 
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I. INTRODUCTION 

Bosonic atoms confined in an optical lattice are a re- 
markably flexible system for exploring many-body physics 
QSSHSBSiiEGil], in which strongly con-e- 
lated physics can be explored, for example, through the su- 
perfluid to Mott-insulator transition lfl2ll . In the superfluid 
regime, a Bose-Einstein condensate exists and experiments 
have explored its properties, such as coherence IH 

anna 

collective modes Il5ll . and transport |0, EH [13] • To date few 
experiments have considered the interplay between the con- 
densate and thermal components that occurs at finite tempera- 
tures lfl3l[l5ll . Indeed, quantitative experimental studies of the 
finite temperature regime have been hampered by the lack of 
an accurate method for performing thermometry in the lattice. 
Recent experimental work has overcome this issue lfl8ll (also 
see IU9IP and finite temperature properties will undoubtedly 
receive increased interest in the near future. 

A unique feature of many-body physics with ultra-cold 
atoms is the opportunity to start from a complete microscopic 
theory and perform ab initio calculations that can be directly 
compared with experiments. In the deep lattice and low tem- 
perature limits, bosonic atoms in an optical lattice provide a 
precise realization of the Bose-Hubbard model [20], originally 
proposed as a toy model for condensed matter physics l2lll . 
However, there is a wide regime of experimental interest in 
which the approximations central to the Bose-Hubbard model 
(nearest neighbor tunneling, local interactions, and neglect of 
excited bands) are not valid. In such regimes it is necessary 
to go beyond the Bose-Hubbard model to furnish an accurate 
description of the physical system. 

Theoretical understanding of the properties of bosons in op- 
tical lattices is still emerging, and accurate modeling is made 
difficult by the combined harmonic lattice potential used in 
experiments, which leads to a complex spectrum, even in 
the absence of interactions I22I l23i l24l 12511 . One approach 
is to use quantum Monte Carlo calculations which, in prin- 
ciple, fully include thermal fluctuations and quantum corre- 
lations. Applications of this approach have mainly been to 
the Bose-Hubbard model l26t 1271 12811 . although recently a 
continuous space algorithm has also been developed for the 
full lattice potential B29I1 . Mean-field methods provide an ap- 



proximate treatment that is much simpler to use and are ap- 
plicable in the superfluid regime where only weak correla- 
tions arise from inter-particle interactions. Extensive stud- 
ies of the harmonically trapped gas have demonstrated that 
the Hartree-Fock-Bogoliubov-Popov (HFBP) mean-field the- 
ory 03011 provides a capable description of thermodynamic 
properties [31], that agrees well with experiments 021 l33ll . 
The development of similar mean-field theories for the lat- 
tice system has been much more limited: HFBP calculations 
have been performed for one-dimensional lattice_systems in 
the continuous HUl and Bose-Hubbard limit J36l[37|l , and 
Duan and coworkers have developed a local density version 
for the three-dimensional Bose-Hubbard model in i38l l39ll . 
To obtain a theory suitable for direct experimental comparison 
over a broad parameter regime, it is necessary to go beyond 
the approach in Refs. [38, 39] to obtain a formalism valid for 
shallow lattices and when excited bands are occupied. 

In this paper we develop a HFBP formalism, based on an 
extended Bose-Hubbard model that includes beyond nearest 
neighbor tunneling, excited band occupation, interactions be- 
tween bands and we discuss an approximate treatment of off- 
site interactions. An important aim of our work is to provide 
a formalism suitable for efficient numerical implementation. 
To achieve this we make use of a local density approximation 
(LDA), that accounts for beyond nearest neighbor tunneling 
and excited bands, and we develop an envelope approximation 
that simplifies the treatment of a general anisotropic harmonic 
confinement to a problem with one independent spatial dimen- 
sion. Combined, the LDA and envelope approximations allow 
us to realize an efficient and practical numerical formulation. 
We show under what conditions it reduces to the simplified 
theory in Refs. f38l l39ll and we numerically investigate the 
features of our formalism. 

In section |TT] we derive the many-body Hamiltonian for 
bosons in an optical lattice with two body interaction, which 
we convert to the extended Bose-Hubbard Hamiltonian. We 
make HFBP mean-field approximations to this in section Hill 
We diagonalize the mean-field Hamiltonian in the LDA, and 
compare our implementation to that of lf38l[39tl in section [TV] 
We derive results on the rich structure of the LDA combined 
harmonic lattice density of states in section|VJ which we com- 
pare to the full diagonalization of the non-interacting Hamil- 
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tonian. In section [V]] we show some important features of 
our numerical implementation and present numerical results 
from our model in section IVIII We compare our predictions 
of thermal properties with results from the full diagonalization 
for the ideal gas and with limited experimental results avail- 
able. We consider the significance of beyond nearest-neighbor 
hopping and excited bands and illustrate the properties of our 
model. In the appendices, we consider the extended Bose- 
Hubbard parameters, including an approximate interpolative 
scheme for off-site interactions. 

II. BOSONS IN OPTICAL LATTICES 

A. Lattice potential and units 

We consider an optical lattice formed by orthogonal stand- 
ing waves, created by two opposing lasers in each direction. 
The laser wavelength Aj (in direction j) is off-resonant with 
respect to an atomic transition. The resulting potential in d 
dimensions, up to an additive constant, is: 



d 

Matt(r) = ^y, sin 2 



(1) 



where Vj is the lattice depth and aj = Xj/2 is the lattice 
spacing in direction j. Most of our results can be generalized 
to the non-separable lattice by adjusting the density of states 
we introduce in section [V] We avoid doing this for notational 
simplicity. 

Except where specifically stated otherwise, our results are 
generally valid for non-cubic lattices and lower-dimensional 
systems. 1 By a cubic lattice, we mean the underlying Bra- 
vais lattice has cubic symmetry (or the equivalent in lower 
dimensions, such as the square case) and that the lattice 
spacings, aj, and depths, Vj, are the same in each ax- 
ial direction. This is the regime of most 3D experiments 

j! d n m m m m m m . 

We will generally present results in recoil units, with the 
unit of length being cij/tt and the unit of energy Er = 



h 2 /8ma 2 where m is the atomic mass and a i 
B. Harmonic-trap potential 



1/d 



Experimentally, atoms are subject to a crossed optical 
dipole [45, 46] potential (due to the focused lasers used to 
make the lattices) and often a magnetic trap also JHH^]. These 
effects are well described by introducing an additional poten- 
tial that is approximately harmonic in form, i.e. 



V tT (r) 



1 



2 2 



(2) 



where u>j is the harmonic trap frequency in direction j. In 
3D experiments, the trap is often spherical or cylindrically 



' However, we do not consider quasi-reduced-dimensional systems, where 
some directions are partially accessible, i.e. kT is of the order of the level 
spacing. 



symmetric (e.g. u> x = u> y ^ lo z ). We consider the general 
anisotropic case in d dimensions. We consider both the lattice 
with Vtr(r) = 0, which we call the 'translationally-invariant 
lattice', and the experimentally relevant combined harmonic 
trap and optical lattice potential, which we call the 'combined 
harmonic lattice'. 

In typical experiments |J3L S3, S^], we find the har- 
monic trapping frequencies to be generally between 2ir x 
18 Hz and 2it x 155 Hz, giving uj/ujr between 0.005 and 0.02 

where u> = Y[j ^>j an d = Er/K is the recoil frequency. 



C. Many-body Hamiltonian 

In this work we consider only bosons, with field operator 
*(r) such that j47ll : 



*(r),£(r') 



= 0, 



#(r),*t( r ') =S(r-r'). (3) 



In the ultra-cold regime, a dilute gas of bosons is described by 
the Hamiltonian 114811 : 



H = / dr* 



fllatt + V tI (r) 



* + | / dr^^M, (4) 



where Hutt = —ti 2 V 2 /2m + Vi att (r), g = 4nh 2 a s /m and 
a s is the s-wave scattering length. 



D. Wannier basis 

We expand the boson field operators in a basis of the Wan- 
nier functions of the non-interacting translationally-invariant 
lattice, Wb(r — Rj), where b is the band index and is the 
lattice site position (see appendix [A), so that we have (as in 
lH): 



%i 



w b (r - Ri 



(5) 



where d b , is the bosonic destruction operator for an atom in 
band b at site i. We note that b and i are discrete d-dimensional 
vectors. For convenience, we shall refer to the ground band as 
b = 0. The Wannier basis is a localized basis for sufficiently 
deep lattices but, for a given lattice depth, there is less local- 
ization for excited bands (see appendix |AV Using a local- 
ized basis significantly simplifies the treatment of interactions 
when off-site interactions are ignored. 

The Wannier states are 'quasi-stationary', since they are not 
eigenstates of H\ a tt, so that there are transitions between the 
different Wannier states in the same band due to the single- 
particle evolution. In particular, the matrix element for hop- 
ping from site R^ to site R; for band b is defined as: 



b.i.i' 



dr wt(v - TLi)Hi att w b (r - R 4 . 



(6) 



There is no inter-band hopping (see dB21 >) with the (non- 
interacting, translationally-invariant lattice) definition of the 
Wannier functions we are using. A change of variables in 
© shows that this formula is dependent on R^ and Ri' only 
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through the difference R; — Rj/. Considering the impor- 
tance of beyond nearest neighbor hopping, we note that the 
ground-band next-nearest-neighbor hopping matrix element is 
as much as 25% of its nearest-neighbor counterpart at Vj = 0, 
but decreases rapidly with increasing Vj, and that beyond 
next-nearest-neighbor hopping is less significant, as shown in 
Fig.[16]in appendix iBl 



E. Extended Bose-Hubbard Hamiltonian 

We now express the Hamiltonian in terms of the operators 
a b i by inserting (0 into dU and we consider the resulting 
terms in this section. 

We assume the trap is slowly varying relative to the lattice 
spacings aj so that: 

J drHr(rK(r-R> 6 Kr-Ri') 



(7) 



where u,; = Vt r (Rj)- I n this work, we will always use the lo- 
cal energy form (0 to represent the harmonic trap. However, 
there are approximations involved in (JTJ which we consider in 
appendixlcl We define the total number operator: 



N = J dr*t( r )^r( r ) =Y f m,i, 

b.i 



(8) 



where h b j = at i a b v Then, expressing the Hamiltonian 
in the grand-canonical distribution to conserve total particle 
number, K = H — fj,N: 



* = E 



b.i 



o E & l,iA 



b 2 ,i2 a b S ,is a b4 M t 1 i 2 i 3 'l 4 

b 1: b 2 ,b 3l b 4 



(9) 



where J7i 1)42 ,i 3 ,i 4 = g J drwt (r-RjJwf (r-R i2 )tu 63 (r— 

b l ,62,03,64 

Ri 3 )wb 4 (r — Ri 4 ). If we restrict to on-site interactions (justi- 
fied in a deep lattice by the Wannier state locality), (O reduces 
to l\ V, l\, where: 



Ki 



E 



E ( Jb ^ ,£L b 



b.i u b,i< 



+ n bti (vi - fi) 



\ E °LA 



b 2 ,i°'b 3 ,i a b i ,iU M.M • 



(10) 



(this interaction term has previously been stated by M50I1 ). We 
retain a smaller set of interaction parameters, i.e.: 



Ubv = 9 / dr \w b (r)w b > (r)|' 



(11) 



which is a good approximation in the typical experimental 
regime, where the interaction parameters are small compared 
to the band-gap energy scale so that we may ignore collisional 



couplings between bands in the many-body state. This ap- 
proximation would need to be revised in the vicinity of a Fes- 
hbach resonance (e.g. see JUl), but this is beyond our scope 
here. 

We derive an approximation scheme for off-site interactions 
in appendix|D] The result is a modification of the interaction 
coefficients. As discussed in appendix ITfl if we use the all- 
site interaction coefficients in our model at Vj = 0, with ap- 
propriate interpretation of the number densities, our model is 
exactly the same as existing no-lattice models. For the non- 
condensate, we find that the effects of off-site interactions are 
negligible for Vj > 5Er. Formulating a consistent theoreti- 
cal description in the shallow lattice limit is fraught for a Wan- 
nier state approach, because these states are delocalized in this 
regime; some work in the shallow lattice has been reported 
1521 . However, our off-site interaction coefficients provide a 
useful interpolation scheme which is accurate in the no-lattice 
case and for moderate to deep lattices. For the condensate, in- 
terference between sites, mediated by the tails of distant Wan- 
nier states, can reduce the interaction coefficient, as discussed 
in appendix |D] All of our work other than appendix [Pluses 
on-site interaction coefficients. 

Other extended Bose-Hubbard work has used various sim- 
plifications of (|9]i: the use of nearest-neighbor hopping and 
nearest-neighbor interactions [49]; the use of ground band 
only, nearest-neighbor hopping and nearest-neighbor interac- 
tions in a homogeneous system 15 311 : the use of ground band 
only and nearest-neighbor interactions OH; and the use of 
nearest-neighbor hopping and on-site interactions in a homo- 
geneous system [55]. 

Limiting to the ground band of a cubic lattice, nearest- 
neighbor hopping (and adding the energy offset Jo,i,i), and 
on-site interactions, the Hamiltonian reduces to the Bose- 
Hubbard model iHHl, which is: 



-J 



E 

(i,i>) 



^ n ,i(vi - Ai) + Y E ™o,i ("cm 



-1), 
(12) 



where restricts the sum to nearest neighbors i and i', 

J = J .i.i>, and U = U 00 . 



m. MEAN-FIELD APPROXIMATION 

A. Mean-field approach: condensate and non-condensate 

We assume that the local number of condensate atoms is 
either macroscopic or zero |T5tL IBUl . so that the field opera- 
tor, \&(r), can be separated into a c-number condensate com- 
ponent (the order parameter), $(r), and a non-condensate 
field operator, ip(r), defined by the usual broken symmetry 
approach, $(r) = (\&(r)\, V>(r) = ^(r) - $(r) so that 



^(r)) = 0. 

The assumption that <£>(r) is a c-number is inaccurate near 
the edges of the condensate, where the local condensate den- 
sity, |$(r)| 2 , is small and just below the critical temperature, 
since fluctuations are important in such regions. 
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We expand the condensate amplitude and the non- 
condensate field-operator in a Wannier basis: $(r) = 
J2i ZiW (r-R t ), ^(r) =E&,i Sb,iW b (r-R^. where we have 
restricted the condensate amplitude expansion to the ground 
band. For an ideal gas this assumption is exact, and with in- 
teractions, the approximation is justified by our assumption 
that interactions are perturbative relative to the band gap en- 
ergy scale. 

From (|5]l and the orthogonality, iA3i , and completeness of 
the Wannier functions (from the completeness of the Bloch 
functions), we get: 



®o i) > 8o,i = a t - Zi,8 b>i = a b i (13) 



for b above the ground band. The operators 5b' a satisfy stan- 
dard bosonic commutation relations. The condensate density 
is: 



|$(r)| 2 = J2 ZiZvvX? ~ R>o(r - R, 



(14) 



allowing for the non-locality of the Wannier states, with con- 
densate number: 

N c = Jdr |<&(r)| 2 = ^|z,| 2 . (15) 



For the non-condensate, we assume that the thermal coher- 
ence length is sufficiently short (long range coherence is ab- 
sorbed by the condensate) that the non-condensate one-body 
density matrix is diagonal in lattice site indices, so that the 
non-condensate density is then given by: 



^ f (r)-0(r)\ = ^2 fi 6,t W r - R; 



(16) 



with hb,i = (sl Jibs )■ The total non-condensate population 



when U/6J > 5.83 at T = fl20|,|6l|]. For typical experimen- 
tal parameters, the transition occurs in 87 Rb when V > 13Er 

(where V = Uj ^/^X but can be v £ lQE R for 23Na S 
(the scattering length of 23 Na is smaller than 87 Rb, and j46Tl 
used a large lattice spacing). The lattice depth for the Mott- 
insulator transition is increased for higher filling factors. 

Making the usual HFBP approximation tf30L l59l l62tl . we 
obtain a quadratic Hamiltonian. Separating this Hamiltonian 
by the number of depletion operators 5\ and 5, appearing and 
by band: 



(18) 



with: 

K A = z 



(1! 

Ki,i = i ^- Jo,i,i'Si',i + Vi - [i 



Uoo \zi\ 2 + 2 Uobfib,i 



(20) 



K 2A i = ^IbAi + (C^ + Ki z ?) > ( 21 ) 

where: 

C-b,i = -~Y]Jb,i,i' Sj' ,j + Vi-fi + 2U Qb \zi\ 2 + 2y^U bb >n b >A, 
i> b' 

(22) 

and Si' a is the shift operator from the site Ri to Rj<, e.g. 

Si f iSb.i — 8b A' . 



C. Gross-Pitaevskii equation 



N = J dr ^ t (r)^(r)^> = ^ n M , (17) B y minimizing the energy functional d (kq^/dz* = 



and we define the b band non-condensate population as Nb 

Ei nb,i- 



using (^Sq A — = 0, we obtain the generalized Gross- 

Pitaevskii equation: 

f- ^ JgAA'Sj'A + Vi - fi + U 00 \zi\ 2 + 2^2Uobn b A z t 



B. HFBP Hamiltonian 



0. 



(23) 



To express the Hamiltonian in terms of the amplitudes Zi, 
and operators, 8b,i, we substitute (fT3l into ( fTOt [59||. How- 
ever, the Hamiltonian still includes up to fourth powers in 
the operators 8b,i- We make a quadratic Hamiltonian sim- 
plification by makin g a mean-field approximation motivated 
by Wick's theorem ll30l l60ll . This is valid in the weakly- 
interacting regime; therefore, our work is not valid in the 
strongly-correlated Mott-insulator case. In a 3D cubic lattice, 
the Mott-insulator transition occurs for the unit-filled system 



We note that if Zi satisfies the generalized Gross-Pitaevskii 
equation, then the terms K\a and K\ . are zero and the next 

contribution comes from K2,b,i- 

When the interaction and trap energy is much more sig- 
nificant than the hopping energy, d23l has the Thomas-Fermi 
solution: 

N 2 = — max ( 0, fi - Vi - 2 ^ U bfib,i ) , (24) 
U °° \ b J 
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where fi is determined by N = J2i \ z i\ + J2b i r ~ l M- 



D. Hartree-Fock 



The Hartree-Fock treatment is obtained by ignoring the 
terms S^zf and S bi z* 2 in #2,6,1 which can then be diag- 
onalized by a single particle transformation, setting 6 bt i — 

j u b,i-j®-b j (where me symbol j indicates a sum over 
modes excluding the condensate). The operators a b ■ are cho- 
sen to satisfy usual bosonic commutation relations: 



a b,j' a b',j' 



= Sbb>Sjj>, 



a b,P a b',f 



= 0, (25) 



and the u&jj modes are an orthonormal basis, i.e. 
E»«M,i u Mj' = satisfying: 



£b,iUb,iJ — EbjUb. 



(26) 



so that J^i K 2 ,b,i = Yl'jEbj&tj&bj- Taking the condensate 
to satisfy the generalized Gross-Pitaevskii equation, we have 
Kq — * -Khf> with: 



-^hf = ^ I - 22 Jo,i,i'Si< 

i \ i' 



^2' E b,3 & l,j\j- 
b,j 



1 ^0° I I 2 I 

+ Vi — fi + — N I Zi 
(27) 



Since the Hamiltonian is diagonal in band b and 
mode j, we can treat the Hartree-Fock modes as non- 
interacting, so that the non-condensate is given by fib,i = 
Ki,j\ 2 n BE {E btj ), where n BE {E) = (e^ - l)" 1 . 



Hi (l M b,i-j| 2 — l u fc-ij'| 2 ) = !• We choose the modes to sat- 
isfy the Bogoliubov-de Gennes equations: 

£b,iUb,i,j + Uobzfvbjj = E b> jU bl i t j, (29) 
£b,iVb,i,j + U Qb z* 2 u b;i ^ = —EbjVb^j. (30) 



The Hamiltonian Kq is diagonal with these solutions |59]: 



+ Vi- fi + — \Zi\ I z t 



Kq — J]] z* - ^2 Jo,i,i'S. 

i \ i' 

b,j \ i ) 



(3D 



and we can treat the quasi-particles as non-interacting which 
leads to: 



n b ,i = ^ (\ u b,i,j\ 2 + \vbj,j\ 2 ) n B E(E b ,j) + \v b ,i,j\" 
j 



(32) 



The references B62L I63H explain that, for a general poten- 
tial, (|29l and ( f30b give quasi-particle functions which are 
orthogonal to the condensate only in a generalized sense, 
J2i z i u b,i,j + z i v b,i,j = 0. To be orthogonal in the sense 
J2i z t u b,i,j — J2i z iV b ,ij = 0, adjustments are required, e.g. 
Eo,jUo,i,j is replaced by ll36ll64ll : 

>■■'• j + Uqo ^ \zi\ (z*u ,i,j - ZiV 0;it j) Zi. (33) 



We do not follow this approach since, in our LDA solution 
below, we approximate by using an orthogonal Bloch form 
for the modes. 



E. Quasi-particle treatment 



IV. LOCAL DENSITY APPROXIMATION 



In general, it is desirable to go beyond the Hartree-Fock 
treatment when the condensate is present, to more fully in- 
clude the effect of the condensate on the excitations of the sys- 
tem (the lattice makes this more important, see section [VTIb . 
To do this, we retain the terms 5 b \zf and 5 2 i z* 2 in the Hamil- 
tonian, which can be diagonalized using a quasi-particle trans- 
formation 15711 : 



S b ,i = ' 



U b,i,3 a b,j + V b,i,j a bJ 



(28) 



The LDA has been extensively used for (non-lattice) har- 
monically trapped Bose gases. The essence of this approxima- 
tion is the replacement — h 2 V 2 /2m — > p 2 /2m in the Hamilto- 
nian with r and p treated as classical variables. The extension 
of this approach to the lattice case is made by the replace- 
ment fii a tt — ► K b (k) where k is the quasi-momentum, b the 
quantized band index and K b (k) the Bloch spectrum. In what 
follows, we present our assumptions in making this replace- 
ment. 



where we refer to the 



o 



as the quasi-particle op- 



erators and the u by ij,v bi ij as the quasi-particle modes. 
We require that d25b holds, as for the Hartree-Fock 

case so that Hi [ u b,i,jUb,i',j ~ v *« « v b-i'.A = and 



b,i,j 



U b,i,j V b,i',j 



V b,i,j U b,i',j 



0. 



The quasi-particle modes are normalized according to 



A. Bloch approximation 

We set j to be the quasi-momentum, k, and make the LDA 
by seeking solutions where u and v have the Bloch form: 

«M',k = e ik ' (Ri ' -Ri Vi,k, Vb.ii M = e ' k ' (Ri '" Ri Vi.k- 

(34) 
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This assumption is exact for the translationally-invariant case, 
and we justify it in general by comparing the non-interacting 
density of states obtained using this approximation to the nu- 
merical diagonalization of the full combined harmonic lattice 
problem in section [VCl 

To make progress, it is useful to consider the Bloch waves, 
V>b,k(r), of ffi att : 



#iatt^i,,k(r) = K b (k)i/j bM (r), 



(35) 



which serves to define the energy, K b (k). We find from (l34l 
and {B3]l that - J2i> Jb,i,i'Ub,i',k = K b (k)u bti>k , so that: 



C. Bogoliubov spectrum 

Making use of the envelope functions from the previous 
section, the Bogoliubov-de Gennes equations, ( f29b and d30l >. 
take the algebraic form: 



£ h (k,r) U 0b z 2 (v)' 
-U ob z* 2 {r) -£ 6 (k,r) 



Ufe(k,r) 
v b (k, r) 



E b (k, r) 



u b (k,r) 
v b (k, r) 
(40) 



Solving the characteristic equation yields: 



C b iU b , 



K b (k)+v i -fX + 2U 0b \z i \ 2 + 



2 5Z U w n b >A u&,i,k- 



(36) 



B. Envelope functions 



We define a function n b (r) which is a proxy with the con- 
tinuous variable r for the number of non-condensate atoms 
per site: n b (R,) = n b: i. Introducing this envelope function 
greatly simplifies our formalism by allowing us to use con- 
tinuous functions to exploit the symmetry of Vt r (r), which is 
broken on short length scales by the lattice. Then, for a suffi- 
ciently small lattice spacing: 

^ / drn b (r) « £n 6 (Ri) =^n 6)i = N b , (37) 

i i 

where a d is the volume of a unit cell of the optical lattice. Sim- 
ilarly, we define the condensate mode envelope z(r), where 
z(Ri) = Zi and n c (r) = \z(r)\ , so that: 

i J drn c (r) » £ |z(R,)| 2 - £ \z t \ 2 = N c . (38) 



We also define the envelope functions u b (k, r) and v b (k, r), 
with u b (k, Rj) = it^i.k and v b (k, Rj) = ify^k, and from 
we have C b: i — > £b(k, r) where: 



A(k,r)=^i,(k)+y tr (r)-/i 

+ 2[/ ob n c (r) + 2^[/^n b ,(r). 



(39) 



Envelope functions represent the discrete functions and do 
not contain the fast Wannier state variation. However, 
apart from exceptional imaging techniques [65], normal op- 
tical imaging techniques would not distinguish density vari- 
ation at the order of one site. If we require the detailed 
spatial density, rather than just site occupation, once we 
have the envelope functions, we can calculate |$(r)| 2 = 
Y, itV z*(R l )z(R 4 ')w5(r - R,)w (r - R*,) from © and 

(^t( r )^( r )\ = J2 b i h b (Ri) Mr - R,)| 2 from ©. 



E b (k,v) = \/^(k,r) - [U ob n c (r)Y. (41) 
From d40l i. choosing the normalization condition | u b (k, r ) | 2 — 



\v b (k, r)| =1 (as in M6CM1 for the no lattice case) we have: 
|2 _ £ 6 (k,r)+£ b (k,r) 



Mk, r)| 
MM)I 



2£ b (k,r) 
2 = £ b (k,r)-£ b (k,r) 
2E b (k,v) 



(42) 
(43) 



Setting v b (k, r) = 0, we find |it b (k, r)| = 1 and E b (k, r) = 
£ b (k, r), yielding the LDA envelope form of the Hartree-Fock 
solution d26l i. 

It has been stated that the Thomas-Fermi approximation is 
necessary to be consistent with the LDA (660 - We use the 
Thomas-Fermi solution for all of our interacting calculations, 
which we restate using the envelope functions, starting from 
to find: 



nJr) = max 



0,M-Mr)-2^E/ 6n 6 (r) 



(44) 



For the non-condensate, using ( l32l and the envelope functions 
we have (BZ is the first Brillouin zone): 



n b (v) = dk 
v 27r/ J bz 



{ 



\u b (k, r)| 2 + |u b (k,r)| 2 n B E[E b (k,r 



(k,r)| 2 }. 
(45) 



From d4TI ). if n c (r) is zero (e.g. above T c or outside the 
Thomas-Fermi radius), we have the Hartree-Fock result. Oth- 
erwise, for the ground band, from d44l >: 



£ (k,r) = K (k) + U 00 n c (r), 



E (k,r) = JK$(k) + 2K (k)U 00 n c (r) 



(46) 
(47) 



which is a useful simplification, and is automatically self- 
consistent with n c (r). 

If we rearrange the equation for the non-condensate enve- 
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lope d45b . we obtain: 
fi ° W = (s)7 Bz dk { 



K (k) + U 00 n c (r) 



E (k, r) 
ifo(k) + £/ 00 n e (r)-£ (k,r) 
2E (k, r) 

2-J J B M E (k,r) C ° th 



n B E[-E'o(k, r) 



/J£fc(k,r) 



-1 
(48) 



If K o (k) is restricted to nearest-neighbor hopping, then this 
result is consistent with that given by Duan and co-workers 
l39tl . We note that they do not make the envelope approxi- 
mation (the discrete LDA sum in their Eqn. (15) should have 
been divided by the number of sites). Additionally, their the- 
ory is restricted to the ground band, and is stated for a cubic 
lattice and a spherical harmonic trap. 



V. DENSITY OF STATES 

The theory we develop relies on detailed knowledge of the 
density of states of the translationally-invariant lattice. 



A. Definition and usage 

By 'density of states', we refer to the per-site density of 
states for the non-interacting, translationally-invariant lattice 
which we define as 16711: 



9b(K) 



1 



(27T)° 



dk8[K - K b (k)], 



(49) 



BZ 



where we take K b (k) from its definition (1351 1. When an inte- 
grand depends on k only through K b (k) we can change vari- 
ables to K = K b (k) since we then have, for any function 

Q b [K b (k),r}: 



dK g b (K)Q b (K,r) 

> 

Applying this to ( |45i >: 

n b {r) = a d dK g b {K) 



(2tt) 



dkQ b [K b (k),r] 



BZ 



(50) 



C b {K,r) 
E b (K,r) 



n BE [E b {K,r)} 



C b (K, r)-E b {K,r) 
2E b (K,r) 



(51) 



We emphasize that this is making no additional ap- 
proximations. Similarly, in the Hartree-Fock ap- 
proach, or above the critical temperature, n b (r) = 
a d dKg b {K) n BE [E b (K, r)]. 

To calculate the density of states, we first need the energy 
dispersion, K b (k), which is easy if the lattice potential is sep- 
arable (the well-studied Mathieu's equation [68, 69l 17011 ). but 
separability is not required. We numerically calculate the den- 
sity of states and show the results in Fig.Q] 




3 4 5 
\K-K (0)]/E R 

FIG. 1 : (Color online) Density of states for the 3D cubic lattice: b = 
(black solid curve), 001 (red dashed curve, the integers specify the 
components b x ,b y ,b z ), Oil (green dashed-dotted curve), 002 (blue 
dotted curve) for (a) V = 5E R and (b) V = WE R 



B. Limiting results for the translationally-invariant lattice 

1. Tight binding 

From ( 1B6I >, the dispersion can be written as a Fourier cosine 
series, with the hopping matrix elements as coefficients: 



K b (k) = -J2 



J 



b i ,3 



2 J2 J Lj cos ( lk i a ^ 

l>0 



(52) 



for a separable lattice, where we define the band b hopping 
between neighbors I sites apart in axial direction j to be • 

(e.g. J l b y = J fc ,ooo,oio and, for the cubic lattice, J = Jqj)- 2 

In the tight-binding limit, beyond nearest-neighbor hop- 
ping is ignored (for the importance of beyond nearest- 
neighbor hopping, see also section I VII BI and appendix [B]). 
In ID, the density of states is then, from d49l ), go(K) = 



V<! 2™| J <l'l 



l-[(K- 



which has infi- 



nite van Hove singularities at the maximum and minimum 
energies of the band, which can also be seen from the zero 



When we use this notation, we are implicitly assuming that the energy 
spectrum is invariant under inversion of quasi-momentum, in view of jBH . 
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> 0.1 0.2 


0.3 0.4 0.5 



[K-Ko(0)]/E R 



[K-K (O)]/Er 



FIG. 2: (Color online) Tight-binding (black solid curve) and actual 
(red dashed curve) density of states for V = 5Er for (a) ID and (b) 
the 2D square lattice 




12 
[K-K (0)]/E R 



0.5 1 

[K-K (0)]/E R 



FIG. 3: (Color online) Tight-binding (black solid curve) and actual 
(red dashed curve) 3D cubic-lattice ground-band density of states for 
(a) V = 2E R and (b) V = 5E R 



derivative in (f52b . In 2D, the square-lattice density of states 3 
has an infinite van Hove singularity at the band center and 
non-zero density at the band edges. The density of states for 
ID and 2D are shown in Fig. [2] In 3D, we compare the tight- 
binding density of states to the actual density of states in Fig. [3] 
for the cubic-lattice ground band. For V > 5Er, the effect of 
beyond nearest-neighbors is much reduced, except for very 
low energies. 



Effective mass 




40 60 
K/Er 



100 



FIG. 4: (Color online) 3D cubic-lattice density of states for V = 
15Er (black solid curve), the free-particle density of states s hifted 
by the minimum energy eig envalue (f ^[K - K (0)]/Er, blue 
dashed curve) and by the spatially averaged energy of the lattice 

(j^(K - \ J2j Vj)/E R , red dashed-dotted curve). 



series, we get the effective mass approximation i\&(k) 



is the effective mass at 
2 & If, 



h 2 k 2 /2m* where m j 
k in direction j, 1/m* = [d 2 K b (k) / dk]] k=ko /ti 
due to the second derivative test, we have m* > for all j and 
assuming that the effective mass approximation applies for all 
K in some region near K™ m (for excited bands and deep lat- 
tices, there is only a small region around ko for which this is 
a good approximation), then for that region of K, from (|49] i: 



9b(K) 



where m* = . 



max (K 



K™ m , 



d/2-l 



T{d/2) (271)^/2 {H 2 /m*) d/2 ' 

f l/d 



(53) 



We note this shows that the van Hove 



singularities at the minimum energy are qualitatively the same 
for the effective-mass assumption as for the tight-binding as- 
sumption: infinite in ID, a finite jump in 2D and an infinite 
derivative in 3D. 

3. High energies 

For high energies, K ^> J2j Vj> me most significant effect 
of the lattice on the density of states is the spatially averaged 
energy of the lattice potential, | . Vj as shown in Fig. [4] 



If, at the minimum energy of a band (K™ m — ^(ko)), C. Limiting results for the combined harmonic lattice 

we have V-Kb(ko) = 0, then from the quadratic Taylor 

In this section, we consider the LDA density of states for 
the combined harmonic trap and optical lattice potential (some 
features of the combined harmonic lattice density of states in 
3 By convolution we can express it as a complete elliptic integral of the first the 1 D tight-binding case, and the 2D case, numerically, are 
kind as g (K)=K i 1- Uk+2 J$ j )/2 J% j ] /4 1 / (2n 2 a 2 Uq J) discussed in lHH). We introduce the LDA density of states for 
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comparison with the full numerical diagonalization as justifi- 
cation of the validity of the LDA approach. 

For the harmonically trapped case, in the non-interacting 
LDA, when we wish to calculate some function, Q[Kj,(ic) + 
Vtr(r)] of the energy, such as the total number of non- 
condensate atoms d37l i and d45l l. we have: 

(2^£ J dr J Bz dkQ[K b (k)+V tr (r)] 
= J dEQ(E) <? L da(£), 



(54) 



from (f50T > where 3lda(^) is given by the convolution: 
9lda(E) ee Y,f dr f Bz dkS i E - J ^( k ) - V ^ 

(55) 



with: 

5tr(Vtr) 



V f dV tI g t r(V tI )g b (E-V tI ), 
b Jo 

(27r) d / 2 



dr5[V tI -V tI (r)] = 



T(d/2) (w 2 ) 



V 



d/2-1 



(56) 



Since the combined density of states, <7ldA (-£■)> nas a rich 
structure, we consider what we expect at various energies. In a 
region where the effective-mass approximation, d53b . applies, 
the contribution to #lda (E) from band b is: 

{d -m^y {E - K ^ )d ~ l > (57) 

where the effective trap frequencies are defined by: 

m 




(58) 



as in JH] and lu* = J] . 0J* 1/d - We therefore expect the initial 
contribution from each band (just after K™ m ) to the combined 
density of states to scale like a harmonically-trapped particle, 
with power d — 1. 

If we assume that the bands are rectangular with width Wb 
and minimum energy K' b nlu , so that gb(K) = l/(W b a ) for 
Wb and gb(K) = otherwise, then: 

2(27r) d / 2 



K mm < R < K mm 



5lda(-B) 



£ 



max (J5 - iq nin , 0) d/ - max (£ - K^ in - W b , 0) rf/2 



W 6 



(2^) d / 2 



r(d/2) (™ 2 af /2 

1 



(59) 

d/2-1 



^g tT [E-Kr n 



2 



(60) 



for E > If™" 1 + Wb, using <[55j and ((56]l. So, we expect 
the eventual contribution of the band to the combined density 



of states (far after K™ m + Wb) to scale like the trap, with 
power d/2 — 1. The high-energy contribution is therefore like 
the density of states for a particle in a harmonic trap with no 
kinetic energy, we call this the 'trap-only' region. 

For energies beyond the effective-mass region, but with 
gmin < E < K min + w ^ me com bi n ed density of states 
depends on the detailed structure of the band g b {K) with an 
approximation given by (1591. 4 

So, the initial contribution from the band is effective-mass 
like and the high-energy contribution from the band is trap- 
only like. We estimate the crossover point between these two 
regimes by equating the single-band contribution from equa- 
tions d57b and (l60l >. In 3D there is no intersection for the first 
excited bands for V > 5Er and, for the ground band: 



F, — K mm 



Wo 
2 



1 



128tH 



,*„2 



h 2 



E, 



jyTain\ z 



(61) 



Using the tight-binding approximations dB7| > and m/ra* w 
7T 2 Jq j/E R j JH] (where E Rij = h 2 /2m\ 2 ) , for the cubic 
lattice and assuming that the cross over is near the middle of 



the band E r 



n 



W /2: 



-E'er 



27 

256tt 5 



W w 0.51 W , (62) 



as shown in Fig. [5] This result has the same scalin g, b ut is 
slightly lower than E a - K Q nin » 0.86 Wo, given in El. 

For high energies, once there have been many bands, we 
consider the assumption that the bands start at the free-particle 
positions, adjusted by the average energy of the lattice (as 
shown in Fig. ©, K™ in = £\ {\Vj + h 2 n 2 b 2 /2ma 2 ). We 
keep the other assumptions leading to d60b and approximate 
the sum in ( f60b by an integral over the region of bands b such 
that < X™ 111 < E, then we recover the density of states for 
a trap with no lattice ( d57| > with m = m*). Evaluating this 
integral in band space, we find: 



Slda(-E) 




i)\{hw) d 



(63) 



so, the eventual contribution of all bands has power d — 1, like 
the density of states of a harmonically-trapped particle. 



D. Comparative results 

We compare the density of states obtained from the full di- 
agonalization of ffi a tt + Vtr(r) (see lf7llD to the LDA den- 
sity of states in Fig. [5] For the low energy LDA results, we 



4 For K^ in < E < K^ in +W b the rectangular assumption implies that the 
contribution to 9lda(^) from band b is proportional to [E — K£ lln ) d ' 2 . 
For 3D, this is a blend between the effective-mass (power d — 1) behavior 
near the start of the band and the trap-only (power d/2 — 1) behavior far 
after the band. For lower dimensions, the rectangular assumption is poor 
from Fig. ff] 
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(E - E^)/E R 



5 10 15 20 

(E - E™ in )/E R 



25 



FIG. 5: (Color online) 3D combined harmonic cubic-lattice density 
of states for (a,b) V = 5E R and (c,d) V = 15E R from the full 
diagonalization (black dotted curve), LDA (red solid curve) and E ct 
(dashed-dotted curve). Shown are [(a),(c)] the ground to the first 
excited bands with the LDA ground band (lower red solid line) for 
reference, [(b), (d)] many bands and the high-energy approximation 
J63b (dashed curve). The LDA is so good that it is obscured by the 
full diagonalization results in all cases. 

also show the contribution from the ground band. We plot the 
product <7LDA(-E)w d , since, for the LDA case, gu^^(E)u d 
is independent of u> from d56l ). For the full diagonalization, 
we can see no dependence of the full density of states mul- 
tiplied by uj 3 for varying u> apart from granularity due to the 
few discrete energies for large to at low energy. The LDA 
results show excellent agreement with the full diagonaliza- 
tion. We note that the approximation ( l63l becomes valid in 
the V = 15E R case only for E > E min + A0E R , beyond the 
region of this plot. The effective-mass region is not visible on 
the plot for V — VSEr due to the scale. 



VI. NUMERICAL IMPLEMENTATION 



A. Translationally-invariant density of states 



£ b (K,f) = K + V %I {f)-iJL+2U 0h n c {f)+2^U w n b ,(r), 

v 

(65) 
(66) 



E, 



■(K,f) = ^£ 2 (K,f)~[U 0b n c (f)} 2 , 
h b (f) = a d l°° dKg b {K) ( l\ n BE [E b (K, f)} 



E b (K,f) 

£ b (K,f)-E b (K,f) 
2E b (K,f) 

We can then calculate the total number using: 

2ir d / 2 



N c = 
N b = 



T(d/2)a d J Q 
2v d / 2 



Y{d/2)a d J 



drr d 1 n c (r), 



drr d 1 n b (r), 



(67) 

(68) 
(69) 



which is now a problem in the two dimensions K and f , and 
is fundamental to our development of an efficient numerical 
algorithm. 

C. Interaction parameters 

We calculate the ID Wannier functions and use their sep- 
arability (from the separability of the Bloch functions) to get 
the interaction coefficients. For the cubic lattice in 3D, the 
densities of the three bands 001, 010 and 100 must be equal, 
i.e. h QQ1 (f) = n 010 (f) = n 1Q0 (f). Thus we can use this sym- 
metry to simplify our calculation of higher bands. For a given 
one of these bands, | of the atomic population is in the same 
band and | is in one of the other first excited bands so that: 

^001,001«00l( f ) + ^ ooi,oionoio (^) + ^ooi,ioo% 00 (r) 

, 0T j > "ooi (r) +n 010 (f) +n 100 (f) 
— (. ^001,001 + ^l / ooi,oioJ 2 , 

(70) 

since J7ooi,oio — C^ooi.ioo- We therefore treat the three ex- 
cited bands together and use (C/ooi.ooi + 2Z7ooi,oio) /3 for 
their self-interaction parameter. 



We find the translationally-invariant energies, K b (k), from 
the non-interacting Bloch solutions to find the density of 
states, by diagonalizing the tri-diagonal (since the lattice po- 
tential is sinusoidal) Hamiltonian, H\ a tt, in momentum space 
rt67ll . We calculate the density of states by binning the ener- 
gies. 

B. Scaled units 



From (|44| | and d45b . n b (r) and n c (r) depend on r only 
through Vtr(r) = \m{w 1 x x 1 + Wyy 2 + u>lz 2 ), so we de- 
fine the scaled co-ordinates x = xcu x /uj 7 y — yu) y /u), z — 
zuj z /u;,f 2 — x 2 + y 2 + z 2 so that V tI (f) = \ 
dxdydz = dxdydz. Our formulae then become: 



n c( r ) =TT max 



0^-V tI (f)-2j2u ob n b (r) 



(64) 



D. Procedure 



We fix the parameters N, Vj , aj , a s , ujj and m throughout 
the entire calculation. For the cubic lattice, we calculate the 
density of states g b {K) and the interaction parameters U b y 
once for each V and use them for any cubic-lattice calculation. 
For the non-cubic lattice, we calculate the density of states and 
interaction parameters for each case. 

We solve d64ll-(lo*7|> self-consistently, finding /.i so that N = 
N c +Ylb^b from (|68T > and j69l . We present our algorithm 
for doing this in Fig. [6] 



Choose it 



Start 
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For each r 



Choose n c (r), set n b (r) = V6 



For each band, b 



Choose n b (r) 



X 



Set: 



C b (K, r) = K + V tI (r) - A* + 2f/ 6n c (f) + 2 £ U w n v (f) 



E b (K, f) = ^(X, f)-[% 6 »ar)] 2 

fi 6 (f) = a d |°° dX fl6 (A') { gff^ "BE[^(g,f)1 



A(-K~,f)-g6(g,f) 
2£ 6 ( AT, f) 



No 



No, next 6 



No, set b = 






r Yes 


Set n c (r) = — — max 
Uoo 


0, f i-V tI (f)-2j2Uobn h (r) 

b 






Yes 



Finish 



FIG. 6: Procedure for LDA calculation 
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We note that, once we have a choice for the chemical poten- 
tial, the calculation is completely local. Therefore, in contrast 
to the Gross-Pi taevskii equation approach of [60], we do not 
check the target for the total number N until the calculations 
at every site are self-consistent. 

For the ground band we use the simplification d48K with 
scaled units and the density of states (this is not shown in 
Fig.®. 

For the translationally-invariant lattice, we use almost the 
same calculation, with V tI (r) set to zero, and use only one 
spatial point, f. However, due to the importance of the low 
energy states in that case, we make the substitution u 4 = K 
and use j dK — > J 4u 3 du so that the integrand isn't diver- 
gent. 



E. Finite-size effect 
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For the non-interacting gas in a combined harmonic lattice, 
we allow for the effect of a positive chemical potential at con- 
densation, equal to the minimum energy /if s = where 
lo* are the effective trapping frequencies, defined in d58b . and 
u* is their arithmetic mean. We limit the domain of the inte- 
gral d67l i to K + Vtr(^) > Mfs> which has a negligible effect 
on results compared to the effect of increasing the chemical 
potential. 

For the interacting gas, it is normal to consider the finite- 
size effect and mean-field interaction shift as independent ad- 
ditive corrections, which we do in [73], but additional work is 
needed to find a consistent way of treating them together. We 
do not consider the finite-size effect due to factors other than 
the positive chemical potential. 

VII. NUMERICAL RESULTS 

In this section we present results demonstrating the ap- 
plication of our mean-field theory to experimentally realistic 
regimes of a Bose gas in a 3D combined harmonic lattice po- 
tential. Our results quantify lattice and interaction effects on 
the thermal properties of the system. We refrain from dis- 
cussing the critical temperature here, which we deal with in 
detail in ifTl . 



FIG. 7: Condensate fraction for a non-interacting combined har- 
monic cubic-lattice in 3D with N = 1000, ui = 0.02uir and 
V — 15Er, comparing full diagonalization (solid curve), LDA with 
A* < Mrs (dashed curve) and LDA with /i < (dashed-dotted curve). 
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FIG. 8: Non-interacting condensate fraction for N = 10 J , w = 
O.OIlor, (a) V = 2E R and (b) V = 5E H . The full diagonaliza- 
tion curve (solid curve) is almost obscured by the all-neighbor result 
(dashed curve) and is appreciably different from the nearest-neighbor 
result (dashed-dotted curve). 

result shows a phase transition (i.e. discontinuous behavior) 
at the critical temperature, whereas the full diagonalization 
shows a more gradual change. 



A. Finite-size effect 

We consider the effect on the non-interacting condensate 
fraction of a non-zero ground-state energy. We plot the con- 
densate fraction for to = 0.02uir and V — 15-Er in Fig. [^re- 
sults at other lattice depths and trap frequencies, are similar, 
except for scaling due to the different critical temperatures). 
We chose a small number of atoms, N = 1000, to accentuate 
the finite-size effect. 

We see that the saturated chemical potential adjustment de- 
scribes the bulk of the finite-size effect well, and the LDA cal- 
culation is in excellent agreement with the full diagonalization 
(by diagonalization of H\ a tt + Vt T (r) to obtain the ideal spec- 
trum which is used solve for the condensate fraction using a 
grand-canonical approach, see [71]). We note that the LDA 



B. Beyond nearest-neighbor hopping 

Here, we consider the effect on the non-interacting conden- 
sate fraction of beyond nearest-neighbor hopping (we use all 
neighbors for our numerical calculations in all other sections). 

We show the condensate fraction for N = 10 5 and lo = 
O.Olwfl in Fig. [8] We see that beyond nearest-neighbor hop- 
ping is significant for V — 2En and much less so for V = 
5Er. For V = 10-Er (not shown), the condensate fractions 
are barely distinguishable on an equivalent plot. The decrease 
in significance of beyond nearest-neighbor hopping with in- 
creasing V/Er, agrees with what we expect from Fig. [3] (see 
also appendixlBl. 
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FIG. 9: Ratio of number of non-condensate atoms in first three (solid 
curve) and beyond first three (dashed curve) excited bands to non- 
condensate atoms in the ground band at the critical temperature for 
the experimental setup of 11211 



C. Excited bands 




15 

V/E R 

FIG. 10: Quantum depletion of 23 Na in a 3D optical lattice. The data 
points with error bars give the experimental quantum depletion. The 
curves give quantum depletion calculated by [46] (solid curve), our 
calculated quantum depletion (dashed curve). 



In this section, we consider the significance of excited 
bands. We do not compare to the full diagonalization, since 
the separation into bands for that calculation is not well de- 
fined. The higher the temperature, the more important excited 
bands are, since they are more thermodynamically accessible. 
We therefore consider the significance of excited bands at the 
critical temperature. It is clear (e.g. see Fig.Q]) that increasing 
the lattice depth decreases the occupation for a given temper- 
ature, and hence the significance, of excited bands. 

We show the number of non-condensate atoms in excited 
bands as a proportion of the non-condensate number in the 
ground band in Fig. [9] The calculations are for 87 Rb using 
HFBP with a s = 5.77 nm and the parameters of lfl2tl with 
an optical lattice wavelength of A = 2a = 852 nm and a 
spherical trap with frequency uo — 2tt x 24 Hz. We used 
their maximum number of atoms, N — 2 x 10 5 . We see 
that excited bands become insignificant for V > "SEr. The 
significance of excited bands at condensation would increase 
for an increased number of particles or a tighter trap, due to 
the increased critical temperature. 

D. Quantum depletion 

The quantum depletion consists of the atoms promoted out 
of condensate due to interactions rather than thermal effects, 
thus leading to a reduction in the condensate fraction at T = 0. 
The number of atoms in the quantum depletion is given by the 
temperature independent part of (05): 

iVQ = 7 ^ F /dr/ Bz dkMk I r)| 2 . (71) 

The quantum depletion is significantly enhanced by increas- 
ing the lattice depth which provides a convenient physical sys- 
tem to explore the crossover from a weakly to a strongly inter- 



acting Bose gas. The experimental measurement of quantum 
depletion in an optical lattice was reported in J46J. In that 
work, atoms were loaded into a lattice, which was linearly 
ramped up to a depth of V ~ 20Er and linearly ramped back 
down. By observing the diffuse background peak of the mo- 
mentum distribution of time-of-flight images during this se- 
quence, the populations of the condensed and non-condensed 
atoms were estimated. The complete ramping procedure led 
to the production of ~ 20% thermal depletion (heating), and 
'Linear interpolation was used to subtract this small heating 
contribution (up to 10% at the maximum lattice depth)' to ob- 
tain the quantum depletion Ir46ll . Their results are presented 
in Fig. [10] We have calculated the zero temperature quan- 
tum depletion to compare with their experimental results. We 
have reproduced their calculations [46] with fixed peak den- 
sity to a level indistinguishable on the plot (solid black curve), 
confirming our microscopic parameters agree with theirs, and 
we found that their results imply N > 10 7 at V = 20E R . 
We used our LDA calculations with fixed total number 5 rather 
than fixed peak density to give improved agreement with ex- 
perimental results with no fitting parameters (dashed curve). 6 
The agreement is improved over the entire range, most no- 
ticeably at higher lattice depths. More precise experimental 
measurements at intermediate lattice depths to better test the- 
ory would be useful. 



5 We have assumed TV = 1.7 X 10 5 atoms, which is mentioned in (4fJ. 
Although the number of atoms throughout is unclear, using their maximum 
number of atoms, TV = 5 X 10 5 , makes only a small change to the results. 

6 We note that our methods are not valid after the Mott-insulator transition. 
Although the n = 1 Mott-insulator transition is at V = WAEji, the 
'measurements were performed at a peak lattice site occupancy number 
~ T Si, and the Mott-insulator transition is at V > 20E R for n > 3, 
which extends our validity regime somewhat. 
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sate plus quantum depletion (dashed curve) and for the Hartree-Fock 
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FIG. 12: (Color online) Spatial densities for the parameters of IT2I1 
at T = 0.8T C , (a) V = 5E R and (b) V = 10E R . Results are for the 
HFBP method: total (black solid curve), condensate (cyan dashed 
curve), quantum depletion (green filled circles); and the Hartree- 
Fock method: total (blue dashed-dotted curve) and condensate (red 
dotted curve). 



E. Effect of quasi-particles 

In addition to the quantum depletion, which was considered 
at zero temperature in section [VII Dl the Bogoliubov quasi- 
particles modify the energy dispersion as in ( 1411 . We compare 
the quantum depletion to the residual Bogoliubov effect in this 
section (using the parameters of 1131 . as discussed in section 
IVIICI ). In Fig. [TT] we show the condensate fraction and the 
condensate plus quantum depletion fraction. At zero temper- 
ature, the only effect of quasi-particles is the quantum deple- 
tion. The methods with and without quasi-particles give the 
same results above the critical temperature and the same crit- 
ical temperature, 7 since equations d66| i and d67l ) are the same 
when there is no condensate. In Fig.QT]we can see the zero 
temperature increase in quantum depletion due to the increase 
in lattice depth (as in Fig. [Tot and we can see that the nature 
of the Bogoliubov quasi-particle spectrum ( f4Tb also increases 
thermal depletion relative to the Hartree-Fock prediction. In 



Fig- 02] we show the total spatial density, and that of the con- 
densate and quantum depletion. The quantum depletion fol- 
lows the condensate density from (|4TT > and d43V A larger lat- 
tice depth increases the effective interaction, decreasing the 
core density and, for the Hartree-Fock case, forces all of the 
thermal depletion away from the condensate region. 



VHI. CONCLUSIONS 

The main purpose of this paper has been the derivation 
of an accurate, computationally tractable theory for describ- 
ing experiments with finite temperature Bose gases in opti- 
cal lattices. Based on an extended Bose-Hubbard model, de- 
rived from the full cold atom Hamiltonian, our theory includes 
the important physical effects needed to describe this system 
over a wide parameter regime. We obtain a mean-field the- 
ory for the system using the Hartree-Fock-Bogoliubov-Popov 
approximation. Through the development of two key tech- 
niques, a local density approximation for the lattice physics 
and an envelope approximation for spatial dependence of the 
mean fields, we realize a formalism for calculation that is ef- 
ficient and accurate. By neglecting the extended features our 
formalism we show that it reduces to a form equivalent to the 
Bose-Hubbard mean-field theory of 13911 . 

We have presented a range of results verifying the accu- 
racy of our theory, and demonstrating the regimes in which 
extended features of our model, over the usual Bose Hubbard 
model, are important. We have also compared to recent exper- 
imental results by the MIT group, and find that our formalism 
provides improved agreement with the experimental data over 
previous calculations 04611 . 

The methods outlined in this paper can be applied to other 
thermodynamic quantities. For example, we have used our 
numerical results to calculate the entropy: 

T =J2 / dr / dKg b (K){(3E b (K,r)n BE [E b {K,v)\ 

b J 

} 



hi 



1 - e 



-PE b (K,r) 



(72) 



and from that the specific heat and then the energy, can be 
obtained. Our formulation is amenable to analytical results as 
we have done in [73]. 

Experimental work in optical lattices is continuing apace 
and, with the recent development of thermometry techniques 
11811 . it is likely that thermodynamics will be measured in the 
near future. For the purposes of developing better understand- 
ing of lattice bosons, and the emergence of beyond mean-field 
effects, it is crucial to have a quantitative and accurate mean- 
field theory for comparison. The theory presented here serves 
this purpose. 
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APPENDIX A: WANNIER FUNCTIONS 

We define the Wannier function for band b, localized at site 
R,- as: 



w b {r - 



R0 = -4= ]T e- ik ' R ^ fc ,k(r), (Al) 



keBZ 



where N s is the number of sites (we let N s — > oo for the 
combined harmonic lattice). We have: 



-0f>,k(r) 



/TV, 



^ e ik ' R ^ 6 (r-R 2 ). (A2) 



For R, on the lattice, X^keBZ e ' k Ri = N s 5n^o, so we have: 
f dru£(r - Ri)wb>(r - R,') = <W<W- (A3) 

For an optical lattice in ID, we show the Wannier function for 
the ground band in Fig. Q~3]and for the first and second excited 
bands in Fig. [14] The harmonic oscillator approximation (the 

eigenstates of Vi a tt(r) ~ 2j=i ( 7rr j / a j) 2 ) overstates the 
peak height at the expense of the tails, and misses the detailed 
structure of the Wannier functions. 

APPENDIX B: HOPPING MATRIX 



Since ifi a tt^6,k(r) = K b (k)ip b ^(r), we have (as in 074J]) 
H iatt w b (r - Ri>) = - J b ,i,i'W b (r - Rj), where hop- 

ping matrix, defined as ([6j: 



6 keBZ 



fc-tai'-ROjjfVk), (B1) 



is the Fourier transform of the energy. In particular, —Jb,i, 
SkeBZ KbQi)/N a is the average energy in the band. So, 



drwl(r — Ri)H laitt w b i (r — Rj<) 
~N, 



- E e - ik ' (R *'- Ri) ^(k), 



(B2) 



keBZ 



so that there is no inter-band hopping and the hopping matrix 
depends only on the difference R^ — R^ . We can invert (IBU 
to write the dispersion relation as a Fourier series: 



K b (k) =-x: Jw • (r " = - e j m, e - ik • r . 

i'=l i=l 

(B3) 

For the ID case, if the spectrum is even in k x then: 

K b Jk x ) = -J bictX - 2^ J L,x cos(lk x a x ). 

l>0 



(B4) 



We demonstrate the Fourier cosine series for the translation- 
ally-invariant lattice spectrum in Fig. Q3] For V = Er, we 
can see that a few terms are needed for the series to approach 
the nearly free-particle dispersion. By V = 5Er, the ground 
band is well described by nearest neighbors. For the first ex- 
cited band, the approach to nearest-neighbor dispersion with 
increasing V/ En is somewhat slower. The width of band b x 
is: 



K h 



= 4 



l>0 



(B5) 



so, for a separable lattice: 



K b (k) = -J2 

j'=i 



J 6°,-,i+ 2 E J ij cos ^' a i) 



;>o 



(B6) 
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FIG. 15: (Color online) Fourier series for the ID translationally- 
invariant lattice spectrum for (a,b) V = Er, [(c), (d)] V = 5Er, 
[(a),(c)] ground band and [(b), (d)] first excited band, using all neigh- 
bors (black solid curve), nearest and next nearest neighbors (red 
dashed curve) and nearest neighbors (blue dashed-dotted curve) 



APPENDIX C: HARMONIC TRAP 

In this work, we will always use the local energy form (0 
to represent the harmonic trap. In this section, we consider an 
exact treatment for the separable case, by defining: 

Vb,v,i,v= [ drV tr (r)w* b (r~R t )w b ,(r-R t ,). (CI) 



1. On-site variation 

Here we consider the accuracy of to the diagonal part of 
Vb,b' There are three components to the integral in (J7j, one 
for each trap direction and the three components are additive. 
Considering, e.g., the x component, we have (using Xi for the 
x component of Rj): 



dr -mUji 2 \wb(r — Rj) 



-vrvJiXf 



da: x 2 \wb(x)\ , 



(C2) 
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FIG. 16: Ratio of beyond nearest-neighbor to nearest-neighbor hop- 
ping. J^j/JIj (solid curve) , J^j/J^j (dashed curve), Jb,j/Jb,j 
(dashed-dotted curve) for (a) ground band and (b) first excited band 



and the width of band b is: 



r/maj; r^min 
K b - K b 



= 4E 



E 

l>0 



J, 



21-1 

bj ,3 



(B7) 



In the tight-binding case where I = 1 dominates, the band- 
width is 4 £.| J*. . 

The ratio of beyond nearest-neighbor to nearest-neighbor 
hopping in shown in Fig.[l6]and we see that the ground-band 
next-nearest-neighbor hopping matrix element is as much as 
25% of its nearest-neighbor counterpart at Vj = 0, but de- 
creases rapidly with increasing Vj. Beyond next-nearest- 
neighbor hopping is less significant. For the first excited band, 
some of the ratios can increase initially. 



since x|-u;b(a;)| is odd and Wb{r) is normalized. For the 
ground band, we can recover (0 by absorbing a constant into 
the chemical potential. For excited bands there is an error due 

to the difference \moj 2 dx x 2 \wb(x)\ 2 — \w (x)\ 2 , 

which is applied to n^, in the Hamiltonian. We plot the con- 
tribution for the first excited band in Fig. [TTJ a). 



2. Off-site contribution 

Now, we consider the case with i ^ i' and b — b'. We 
note again that the components of the trap contributing to the 
integral in the three directions are additive. We only get a 
potential error in the x component if components in the other 
directions of i and i' are equal. Then, for Xi ^ Xy : 



J dr ^mw 2 x x 2 w* b (r - Rj)u> b (r - Rj 



= -mu , / dxx 2 wl(x)w b (x - (X v - Xi)), (C3) 

since Wb (x — Xi) and Wb (x — Xy ) are orthogonal and (x — 
Xi)wb(x — Xy) is even about (Xi +Xj')/2 as Wb(x) is either 
even or odd. In Fig. fTTl b) we plot this contribution for nearest 
neighbors as a function of V. 



3. Inter-band contribution 

Now we consider the case with b ^ b' and i = i'. To allow 
for this contribution, it would be necessary to include matrix 
elements between bands in the Hamiltonian. 

To quantify the error, we consider the additive component 
in the x direction. There is only a contribution if the other 



17 



0.5 
0.4 
0.3 
0.2 
0.1 


0.4 
0.3 
0.2 
0.1 


0.2 

0.1 





(a) ■ 




\ 




\ 


(b) ■ 


\ 




\ 




s 













10 



15 



V/E R 



FIG. 17: Error due to (a) on-site variation of the trap, I\ — 

f^° oa dx(x/a x ) 2 [\wi(x)\ 2 — |wo(a;)| 2 ] (this form is chosen so that 

the error is in units of ^muj 2 a x ), (b) contribution from adjacent 
sites, I2 = j da; (x/a x ) 2 wl(x)wt(x — a x )\, ground band 
(solid curve), first excited band (dashed curve), (c) Wannier func- 
tion overlap between the ground and first excited bands, I3 — 

dx (x/a x )wo(x)w 1 (x). 



components of band b and b' are equal. Then, with b x being 
the x component of b and b x ^ b' x : 

J dr ^mujlx 2 wl(r - Ri)w b >(r - R. t ) 

dx (a; + Xi) 11}% (x)wb> (x). 



(C4) 



Considering, e.g. b x = (the ground band), and b' x = 1 
(the first excited band) Wq(x)wi(x) is odd so the above be- 
comes mL0 x Xi Axxwq{x)w\(x). In Fig. HTl c). we plot 
this contribution as a function of V . 



APPENDIX D: INTERACTION COEFFICIENTS 

1. Beyond the on-site interaction approximation 

Here we derive approximate results for interactions extend- 
ing to all sites. To do this, we make the HFBP mean-field ap- 
proximations, as discussed in section|III] but starting from the 
more general extended Bose-Hubbard Hamiltonian ©. As 
in the on-site case, we ignore collisional couplings between 
bands in the many body-state. For the non-condensate, we 



also ignore collisional coupling that relies on coherences be- 
tween sites (i.e. requiring two indices at two sites) in the 
many-body state, to find: 



*1 >»2 >*3 >M 
b l > b 2 > b 3 ' b 4 

* h h' h h' 



b l ' b 2 > b 3 > b 4 



(Dl) 



.,6,6' 



We assume that the density varies sufficiently slowly that 
nb,i ~ ^6,3 for Sltes Rj near Ri' In the following, we will 
sum over all sites, by assuming that where the approximation 
Ub,i ~ nb,j is poor, due to the sites being far apart, these terms 
will be suppressed by the negligible Wannier function overlap. 
Then we have: 



*1 > 4 3 >*4 
b l > fa 2 > fa 3 > fa 4 



b l .&2i b 3= b 4 



s» Ag ^2 ^6,i^6,i%,i J dr \ w b( r ) w b> (r - Ri')l 

i,b.b' i' 

= 4 Y KJb,inb>,iUlv, 

i,6,6' 



(D2) 



which is the same as in (1211 with U' w substituted for Ubv 
where: 

U W=9^2 J drltflfc^tflyfr-RiOl 2 . (D3) 

For the coherent condensate, we assume that w Zj, for sites 
TLj near R^. As above, we assume that contributions between 
sites far apart are suppressed by the negligible Wannier func- 
tion overlap. Assuming that the phase factors are chosen so 
that wo (r) is real, we have, for site R^ : 

= 9 Y z *i z i2 Z is z i4 / dr n^o(r-Ri 3 .) 

12,13,14 ' 3=1 

4 f 4 



3=1 



= g\z il \ 4 J dru. (r) [0vTVo,o(r)] 3 - |z 4l | 4 [/^ , (D4) 

where X)i w &( r ~ -^») = V / ^YsV'&,o( r ) is the Bloch func- 
tion normalized over a single site, from (1A21 >. and ?/>o,o(fi) oc 
ceo(rj7r/aj, g) (the Mathieu function) is real and periodic on 
the lattice. The result takes the same form as above with Uq 
substituted for [/no where: 



Udo =g / drwo(r) V^ ,o(r) 



(D5) 



Similar arguments could be used for the terms involving inter- 
actions between the condensate and the non-condensate. The 
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FIG. 18: Interaction coefficients in 3D. (a) Ground band. On site 
(dotted curve). All sites: non-condensate U' iD3\ (solid curve), con- 
densate U" JD5b (dashed curve). No lattice limit: all sites JD6b (x), 
on site (Q77J (o). (b) Excited-band. On-site: 000,001 (these inte- 
gers specify the components b x ,b y , b z of each band) (solid), 001,001 
(dashed curve), 010,001 (dashed-dotted curve); corresponding all 
sites (dotted curve). 



above results are appropriate for the pure thermal gas, e.g. for 
finding the critical temperature from above, and for the pure 
condensate at zero temperature. To quantify the effect of off- 
site interactions on the thermal depletion, terms for interac- 
tions between the condensate and the non-condensate would 
be needed. 



2. No lattice limit 

When there is no lattice, the Hamiltonian fl35l l gives us 
Kh(k) = h 2 k 2 /2m and the Bloch states are plane waves. 
Using these to evaluate the Wannier functions from ( IA1I ). 
and then the all-sites interaction coefficients: C/q easily from 
dD5l l, and U' w from flD3t by splitting the sum into axial com- 



ponents and recognizing the Riemann zeta sums to get: 

^00 = ^00 = ^000,001 = ^001,001 = ^010,001 = 



(D6) 



So that, if we use all-site interaction coefficients and also 
treat n c (r) and n b (r) as the condensate and non-condensate 
densities (rather than as envelope functions, with densities de- 
fined by (TBI and ( fToT i. although the total condensate and non- 
condensate numbers do not depend on this distinction, from 
([T5l > and ([T7| >) then all of our LDA equations in section [TV] 
would be the same as we would get from a no lattice calcu- 
lation [60], in spite of our expansion of the field operators in 
a Wannier basis. When only on-site interactions are included 
there is a shortfall, using (fTTT i: 



U = -Qj ( — ) , £/ooo,oon — 



a 3 27' 

u 9 2 u 9 25 

<-^00n,00n — ^3 g ' U 0n0,00n — —3 ' 



(D7) 



of, for example, 1 - (2/3) 3 = 70% for the 3D ground-band 
coefficient. For reference in Fig. [18] a/E^a s = 8a 3 /gir. 
3. Comparison 

The 3D ground-band interaction coefficients are shown in 
Fig. [T8l a). Both all-sites interaction coefficients, Uq and Uq , 
include their corresponding on-site component, f/oo, in their 
sums, (ID31 > and (D4i . For the non-condensate interaction co- 
efficient, all other terms in the sum are positive (since we have 
excluded interference), so that off-site interactions always in- 
crease the interaction coefficient (relative to U). 

The 3D excited-band interaction coefficients are shown in 
Fig. [T8T b). The results all tend to the expected limits at V = 0. 
The gap between all-site and on-site interaction coefficients is 
maintained for higher V/Eji than for the ground-band, since 
the excited-band Wannier functions are less localized. 



APPENDIX E: DIAGONALIZATION OF THE QUADRATIC HAMILTONIAN 



This appendix gives a derivation of the quadratic Hamiltonian HSi . and a proof that the Bogoliubov-de Gennes equations 
reduce the quadratic Hamiltonian to diagonal form (l3TT >. 



1. Quadratic Hamiltonian 



We begin with the extended Bose-Hubbard Hamiltonian: 



E 



f>l,i>2,f>3>&4 



b l ,&2! b 3>&4 



(El) 
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We make the substitutions Zi = (^a Q { j , 8qj = a Q i — for the ground band, and 8 b) i = a b i above the ground band (with the 
operators 8 b .i satisfying standard bosonic commutation relations) into the interaction term of ( IEU to obtain: 

*rt/00 + > ]lz!Sb.i + Zi6t„) \Zi\ 2 Ut.i.<,< 



61,62,63,64 

+ E Q< 2 <M6<,i + ~*?$A< + 2 N 2 *JA') ^;-*.'-«, + E (*? ^ A A,< + ^AA*) ^ t M.M t 

6,6' ^ / ' • • 61,62,63 1,2,3 

+ ^ E ^.AAA,^ M.M , (E2) 

Z z * ' ' °1,°2- D 3> b 4 

61,62,63,64 

where we have assumed phase factors are chosen so that the Wannier functions are real, so that the order of subscripts in 
U i.i.i.i is unimportant. 

fe l> fa 2, fa 3> fa 4 1 — -1 

We make a quadratic Hamiltonian simplification by making a mean-field approximation motivated by Wick's theorem 1 30 ] . 
For the fourth order terms, we find: 

\ E ^AAA,^. W, -^y^nb'JtJb^Uw, (E3) 

Z * — ' 01,02,^3,64 * — » 

61,62,63,64 6,6' 

where Uw = g J dr \w b {r)wb' (r)\ 2 and we have used a Popov approximation to eliminate the terms /j| ^ ^ and (^8b,i8 b ^, 

and we neglect pairs with different band indices, since we ignore collisional couplings between bands in the many-body state. 
Similarly, we simplify the third order terms by analogy with Wick's theorem M62I1 to find: 

E ih 2 ,ih 3 ,iU « 2z*8 ,i y^n bti U b, (E4) 

' ' ' 0,bi,6 2 ,D3 ' ■ 

61,62,63 6 

and the adjoint of this equation. We set the linear terms iz*5b,i + z i$\ %j \ z i\ to zero f° r ^ 7^ and the quadratic terms 



|zi| 2 5| j^6'.i, z* 2 8b,i8b',i and z 2 2 <5| ^ 5^, 4 to zero for b 7^ b' by the same assumption that interactions are perturbative relative to 
the band-gap energy scale. Our interaction term becomes: 

\ E "L* a L%,z%,* U b *•*•*•* b ~ U + + ^o,i) ki| 2 ^00 

61,62,63,64 J ■ 1 \ / 

+ (\ z ^^h + \ z f 8 Z + 2 l^l 2 + 2z*n bii 8 ,i + 2z i n b , i 8U U ob + 2^%,i^/ M ^6', (E5) 

6 ^ ' 6,6' 

which gives: 

= E (^M + #M + k h + E *3,mJ ' (E6) 



with: 



where: 



K .i = z* ^- ^ Jo,i,i>Si>,i + Vi- fi+ N^j z *> ( E7 ) 

K lti = 5~l ti f- ^ Jo,i t i'Si^i + Vi - [i + Uqq \zi\ 2 + 2 E ^06"6,i^ Zi, (E8) 

& w ^ siAA, + u ~y {}fA + Ki*f) . ceo 

A>,« = ^ E J b,iA>Si> yi +Vi-fi + 2U ob \zi\ + 2 Uwfiy ii , (E10) 



6' 



and Sj^j is the shift operator from the site R; to R;<, e.g. t i8 b: i — 8 bt i< . 
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2. Quasi-particle treatment 



The quasi-particle transformation 



(Ell) 



(with the operators a b satisfying standard bosonic commutation relations) brings Ki,b,i mt0 tne f° rm: 



j,k 



a b,j a b,k 



a b,j a b,k 



l b,i,j£b,iUb,i,k + ~Y ( z i u b,i,j v b,i,k + z * 2 vl^jU bt i^ k ) 

A * . Uob (2 * . *2 * \ 

v b,i,j!~b,iV b i k H — \^Z i V b jjU bik + Z i Ub,i,jVb,i,k) 

r i ^06 / 2 *2 \ 
Vb,i,jLb,iUb,i,k + — { z i Vb,i,jVb,i,k + Z i Ub,i,jUb,i,k ) 

* A * i Uftb / 2- * * , *2 * * \ 

u b,i,j J -b,i v b,i,k "I ^~ \ z i u b,i,j u b,i,k "+" z i v b,i,j v b,i,k) 



(E12) 



To calculate the tunneling term, we first consider a property of the shift operator, SV^. Since Jb,,,^ = J^/ we have: 

^ l ^i Jb,i,i' 8j' jUi = ^ ^ (^Jb,i,i' Si> iXi*j l/i, (E13) 



by interchanging the roles of the dummy variables. 8 From (lElOK since the diagonal terms in £ are real, we therefore have 9 
J2 % x*Cb,iVi = J2i \ £b,iXi) Vi so that: 



^xlCb.iVi 



^2 \ £b,i%i) Vi +^2x*£bal 



(E14) 



and: 



E^,M = ^E' [i E bd + E b ,k) (&l 



At 



i,3,k 



\j a b,k u b,i,j u b,i,k ~ 0i bij a h k v b ,ijv b<iik 

,., - • lb,k) f "(,.,•"/../, ' '■•'■./ 
We choose the modes to satisfy the Bogoliubov-de Gennes equations: 



-{E b ,j - Eb,h) (& b .i & b,k v ^i,J u hi,k ~ &bj & b,k u b,i,j v h,k 



£b,iUb,i,j + Uo b zfv byi j = Eb t jUb,i,j, 
£b,iV b ,ij + Uo b Z* 2 U by ij — —EbjVb,i,j- 



(E15) 



(E16) 
(E17) 



The second term in ( |E15t is directly zero for j = k and from v b: i. k x ( |E16t + Ub,i,k x dE171 > and applying £ to the left: 

(E bJ + E b , k )(u btiJ v b j ik - Vb,i,jUb,i,k) = 0, (E18) 



so, for j ^ k we have: Vb,i,jUb,i,k — Ub,i,jVb,i,k> taking E bl k to be non-negative II58I1 . Therefore, the sum of each pair of 
opposite off-diagonal elements of the coefficients of a b j& b k is zero. The same argument works for the off-diagonal coefficients 
of a\ jCt b k using the complex conjugate. 



This result continues to apply if we exclude, e.g. beyond nearest or beyond 
next-nearest neighbors by symmetrically setting J b i i i = for hopping 
terms not required. 

9 This result shows that C is Hermitian under the inner product (x\y) = 
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The first term of ( IE15b becomes: 

(Eb t j+Eb,k) al o a b k (ub ti jUb,i,k-Vb,i,j v b,i,k)~^k \Vb,i,j\' 



(E19) 



where we have exchanged the dummy variables j and k for the a b k terms. From ul { k x dE16b + v b i k x jE17b and applying 
t to the left: 



(E b ,j - E b ^ k ){u b ^ tj ul hk - v b ,i,jvt i k ) = 0, 



(E20) 



so taking the complex conjugate for j ' ^ k we have i it^i.fc = v bi .Ub^fc, eliminating the off-diagonal terms, and using 



e/(k 



\ v b,i,j\ 



= 1 for the diagonal terms, the Hamiltonian is reduced to the diagonal form: 

U 00 



KQ = Yl z i [~Y1 J o,M"SV,* + Vi - M+ ~y \zi\ 2 J z t + ' E b ,j I a\ tj a bd - ^ \v b ,i,jf j 

i \ i' / b,j \ i / 



(E21) 



r 
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